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A large part of the interest in magnets with frustrated antiferromagnetic interactions comes from 
the many new phases found in applied magnetic field. In this Article, we explore some of the new 
phases which arise in a model with frustrated ferromagnetic interactions, the J 1 -J 2 -J 3 Heisenberg 
model on a square lattice. Using a combination of classical Monte-Carlo simulation and spin- 
wave theory, we uncover behaviour reminiscent of some widely-studied frustrated antiferromagnets, 
but with a number of new twists. We first demonstrate that, for a suitable choice of parameters, 
the phase diagram as a function of magnetic field and temperature is nearly identical to that of 
the Heisenberg antiferromagnet on a triangular lattice, including the celebrated 1/3-magnetisation 
plateau. We then examine how this phase diagram changes when the model is tuned to a point where 
the classical ground-state is highly degenerate. In this case, two new phases emerge; a classical, 
finite-temperature spin-liquid, characterised by a “ring” in the spin structure-factor tS(q); and a 
vortex crystal, a multiple-Q state with finite magnetisation, which can be viewed as an ordered 
lattice of magnetic vortices. All of these new phases persist for a wide range of magnetic field. We 
discuss the relationship between these results and published studies of frustrated antiferromagnets, 
together with some of the materials where these new phases might be observed in experiment. 

PACS numbers: 75.10.Hk, 67.80.kb, 74.25.Uv 


I. INTRODUCTION 

Much attention has been devoted to the question of 
whether a frustrated magnet orders or not^. Even in the 
cases where such systems do order, the results are often 
surprising and unconventional. Novel phases of matter 
have often been uncovered by applying an external mag¬ 
netic field to frustrated magnet systems^’^. Examples 
of these phases range from quasi-classical magnetisation 
plateaux and magnetic supersolid phases^"^, magnetic 
analogues of liquid crystals^”^^, and crystals formed of 
skyrmions, which have already been observed in experi- 
ment^^’^^. 

This ongoing interest in novel magnetic phases is un¬ 
derwritten by a steady stream of new materials. Among 
these, materials with frustrated ferromagnetic interac¬ 
tions, i.e. with ferromagnetic largest interactions but 
with a non-ferromagnetic ground state driven by compet¬ 
ing interactions, have a particularly rich phenomenology. 
Well-known examples include spin ice, a celebrated ex¬ 
ample of a three-dimensional classical spin liquid and 
the solid phases of ^He, whose nuclear magnetism contin¬ 
ues to push the limits of our understanding of quantum 
spins^^’^^. More recent discoveries include new families of 
layered, square-lattice vanadatesand cuprates^^“^^, 
with properties that encompass ordered ground states 
selected by fluctuations^^, exotic singlet phases^^, helical 


order^^, and surprisingly, a 1/3-magnetisation plateau^^. 

In this Article we explore some of the novel phases 
which arise in a simple example of a frustrated ferromag¬ 
net — a Heisenberg model on a square lattice, in which 
ferromagnetic 1^^-neighbour interactions compete with 
antiferromagnetic 2^^- and 3^^^- neighbour exchange. 
Using a combination of classical Monte Carlo simulation 
and spin-wave theory, we establish the phase diagram of 
this model as a function of temperature and magnetic 
field, for two different sets of exchange parameters. In 
the process, we uncover a number of phases not usually 
associated with square-lattice magnets. These include 
a collinear 1/3-magnetisation plateau, a spin liquid with 
finite magnetisation, and a crystal composed of magnetic 
vortices. Illustrations of these phases, and of the phase 
diagram for one of the parameter sets, are shown in Eig. 1 
and Eig. 2, respectively. 

The model we consider is a natural generalisation of 
the “J 1 -J 2 ” model, which describes competing Heisen¬ 
berg exchange interactions on the I^^- and 2’^^-neighbour 
bonds of a square lattice. The J 1 -J 2 model has a long and 
interesting history, and remains one of the paradigmatic 
examples of a frustrated magnet^^. Much of the original 
interest in this model was driven by the possibility that 
competing antiferromagnetic exchange could stabilise a 
quantum spin-liquid, a question which continues to in¬ 
spire new research^"^"^^. However the J 1 -J 2 model is also 
significant as a working model of magnetism in iron pnic- 
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FIG. 1. (Color online). Novel phases found in a square-lattice frustrated ferromagnet in applied magnetic field, (a) Collinear 
1/3-magnetisation plateau, with stripe-like three-sublattice order, (b) Spin liquid with short-range helicoidal correlations, 
characterised by a “ring” in the spin structure factor <S(q), shown here for /i = 0. (c) Crystal of magnetic vortices, with finite 
magnetisation. The colour map is a representation of the phase of the spin texture in the S^-S^ plane. The model studied is 
the J1-J2-J3 Heisenberg model on a square lattice, [Eq. (1)]. 


tides^^“^^, and as one of the simplest possible examples 
of a frustrated ferromagnet. In particular, where ferro¬ 
magnetic 1^^-neighbour exchange competes with antifer¬ 
romagnetic 2^^-neighbour exchange, the spin-1/2 J 1 -J 2 
model can support spin-nematic order^’^^“^^. 

Extending frustrated spin models to include further- 
neighbour interactions typically results in classical 
ground states with incommensurate, helical order. For 
large values of spin, the leading effect of quantum and/or 
thermal fluctuations is a small correction to the pitch of 
the helix, leading to small, quantitative, changes in phase 
boundaries^^”^^ . In the case of the J 1 -J 2 -J 3 Heisen¬ 
berg model on a square lattice, these expectations are 
borne out by linear spin-wave theory^^”^"^. However ex¬ 
act diagonalisation calculations for S = 1/2 suggest a 
richer phase diagram, including a number of novel ground 
states^^’^^. And even at a classical level, other models, 
such as the J 1 -J 2 -J 3 Heisenberg model on a honeycomb 
lattice, can support a much richer behaviour"^^’^^. 

One of the possible alternatives to simple helical order 
are multiple-Q states, composed of a coherent superpo¬ 
sition of states with different ordering vectors. In gen¬ 
eral, these states are not compatible with the fixed spin 
length constraint |S| = 1, and therefore do not belong to 
ground-state manifolds. However, there are a number of 
recent interesting cases where these phases are stabilised 
by the interplay between frustration and thermal fluc¬ 
tuations, such as the classical pyrochlore antiferromag- 
net^^ or the classical extended triangular-lattice antifer- 
romagnet^^. Interestingly, multiple-Q states also arise 
in the description of spin textures composed of crystals 
of topological defects. These are the celebrated case of 
skyrmions lattices in chiral magnets such as MnSi^^’^^, 
magnetic-vortex lattices in a generic class of Mott insu¬ 
lators^^, or skyrmion crystals in the antiferromagnetic 
triangular lattice^^. Magnetic-vortex lattices in insulat¬ 
ing systems seem to be especially uncommon outside the 
realm of superconducting systems^^’^^. 


The picture which emerges from the square-lattice 
frustrated ferromagnet studied in this Article is inter¬ 
esting for a number of reasons. Just as in frustrated 
antiferromagnets, the interplay between ground state de¬ 
generacy and fluctuations leads to a rich behaviour in 
applied magnetic field — cf. Fig. 2. However, unlike 
those antiferromagnets where the frustration comes from 
the geometry of the lattice, the model under study here 
can be tuned to exhibit a wide range of different mag¬ 
netic phenomena, many of them unexpected on a square 
lattice. This is particularly true when parameters are 
chosen so to place the model at the border between com¬ 
peting forms of order. 

The remainder of this Article is structured as follows: 

In Sec. II we introduce the model and the methods used 
to study it. On the basis of known results for the clas¬ 
sical ground-state we single out two parameter sets; one 
which corresponds, at a mean-field level, to a triangular- 
lattice antiferromagnet; and another for which the clas¬ 
sical ground-state is highly degenerate. 

In Sec. Ill we present Monte Carlo simulation results 
for the first of these parameter sets. We demonstrate that 
the phase diagram as a function of magnetic fieid and 
temperature is aimost identicai to that of the Heisenberg 
antiferromagnet on a trianguiar iattice, and inciudes its 
ceiebrated 1/3-magnetisation piateau, now transiated to 
three-subiattice order on a square iattice. These resuits 
are summarised in Fig. 2. 

In Sec. IV we present Monte Cario simuiation resuits 
for the second of these parameter sets. We estabiish the 
phase diagram of a finite-size ciuster as a function of 
magnetic fieid and temperature, and present evidence 
for both a finite-temperature ciassicai spin-iiquid and, 
at iow temperature and iow vaiues of magnetic fieid, a 
muitipie-Q state with the character of a vortex crystai. 
These resuits are summarised in Fig. 8. 

In Sec. V we discuss how these resuits reiate to pub- 
iished work on frustrated antiferromagnets, how quan- 
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FIG. 2. (Color online). Phase diagram of a square-lattice 
frustrated ferromagnet in applied magnetic field, exhibiting 
the same phases as found in the Heisenberg antiferromagnet 
on a triangular lattice [57] — a coplanar “Y-state”, inter¬ 
polating to “120-degree” order at vanishing magnetisation; 
a collinear 1/3-magnetisation plateau [cf. Fig. 1(a)]; and a 
coplanar 2:1 canted state. Phase boundaries are taken from 
classical Monte Carlo simulations of [Eq. (1)], for the 

parameter set Ja [Eq. (18)], and scaled to the thermodynamic 
limit, as described in Section IIL Temperature and magnetic 
field measured in units of Ja [Eq. (13)]. 


turn effects might enter into the problem, and where these 
novel phases might be realised in experiment. 

Finally, in Sec. VI, we conclude with a brief summary 
of the results and remaining open questions. 

Technical details of associated spin-wave calculations, 
used to confirm the results of Monte Carlo simulations, 
are discussed in a small number of Appendices. 

Appendix A develops the mathematical formalism 
needed to carry out both a classical low-temperature ex¬ 
pansion, and a linear spin-wave expansion, about the 
classical ground states found in this model. 

Appendix B applies this analysis to the ground state 
manifold for the second parameter set, and shows that 
thermal fluctuations select a spiral state at low tempera¬ 
tures and zero field, as suggested from Monte Carlo sim¬ 
ulations. 

Appendix C analyses how, for the second parameter 
set, a canted Y state eventually prevails over a conical 
version of the spiral state when a magnetic field is ap¬ 
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FIG. 3. (Color online). Exchange interactions up to 
3”^-neighbour on the square lattice. In this Article we con¬ 
sider a Heisenberg model [Eq. (1)] where ferromagnetic 

1^^-neighbour exchange Ji competes with antiferromagnetic 
2'^^- and 3”^-neighbour exchanges J 2 and J 3 . 

plied. The vortex crystal is not favoured at very low 
temperatures, again in agreement with the Monte Carlo 
results. 

Appendix D discusses, from linear spin-wave calcu¬ 
lations of the quantum corrections to the ground-state 
energy, how the classical phase diagram for the second 
parameter set changes in the quantum model at low mag¬ 
netic fields and temperatures. 

II. MODEL, METHOD AND ORDER 
PARAMETERS 

A. The J 1 —J 2 —J 3 Heisenberg model on a square 
lattice 

In this Article, we consider one of the simplest possible 
prototypes for a frustrated ferromagnet, the Heisenberg 
model on a square lattice, with competing 2'^^- and 3^^^- 
neighbour exchange 

= Ji ^ Si • S,- + J 2 ^ Si • s,- 

(u>l (u>2 

+j35]Si-s,-/i^5y (1) 

(U)3 * 

Here is a classical spin with = 1, the sum 

runs over the n^^-neighbour bonds of a square lattice, as 
illustrated in Fig. 3, and h is an applied magnetic field. 
We restrict ourselves to the case where 1 ^^-neighbour ex¬ 
change Ji is ferromagnetic, while further-neighbour ex¬ 
change J 2 and J 3 are antiferromagnetic, i.e. 

Ji < 0 , J 2 > 0 , J 3 > 0 . ( 2 ) 

The exchange integrals Ji, J 2 , J 3 define a 3-dimensional 
parameter space, and for many purposes, it is convenient 
to represent them as a vector 

J = J 2 , J 3 ). (3) 

The classical ground-state phase diagram of 
together with its spin-wave excitations. 
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was studied in a series of papers by Rastelli et 
and Chubukov^^. In the absence of magnetic field, 
allowing for the leading effect of fluctuations, all ground 
states are found to be coplanar spirals characterised by 
a wave vector 

Q = {Qx^ Qy) • ( 4 ) 

This wave vector can be determined by minimising the 
Fourier transform of the interactions 

J(q) = 2 Ji(cosga, + cos qy) + 4J2 cosg^, cos Qy 

+2J3(cos + cos 2g^) . (5) 

For ferromagnetic Ji [cf. Eq. (2)], there are four distinct 
cases : 

1. a uniform ferromagnetic (FM) phase with 

Q™ = (0,0). (6) 

2. a two-sublattice collinear antiferromagnetic (CAF) 
phase with 

qCAF ^ (^^0) or (0,7r). (7) 




FIG. 4. (Color online). Classical ground-state phase 
diagram of the frustrated Heisenberg ferromagnet on on a 
square lattice, [Eq. (1)], as a function of competing 

antiferromagnetic exchange J 2 and J 3 , following Ref. 42. 
Ground states in the absence of magnetic field comprise 
a collinear ferromagnet (FM), a collinear antiferromagnet 
(CAF), a one-dimensional coplanar spiral (ID spiral) and 
a two-dimensional coplanar spiral (2D spiral), as defined in 
Eqs. (6-9). The two parameter sets studied in this Article, 
Ja [Eq. (18)] and Jb [Eq. (22)], are labeled with black dots. 
A blue dashed line indicates parameters for which the ground 
state is a ID spiral with three-sublattice order [Eq. (10)]. 


3. a family of one-dimensional (ID) spirals with 
= (eiD,0) or (0,(3id) , 


cosQid = - 


Jl + 2J2 

4 J 3 


Jl 


2 J 2 + 4 J 3 


( 8 ) 


4. a family of two-dimensional (2D) spirals with 

= (Q 2 D, Q 2 d) , 
cos(52D = 


(9) 


These ground states, and the range of parameters for 
which they occur, are illustrated in Eig. 4. 


B. Three—sublattice order on the square lattice, 
and its connection with a triangular—lattice 
ant iferromagnet 


In general, the one-dimensional spiral ground states of 
incommensurate, with a wave vector 
[Eq. (8)] which interpolates smoothly from one phase 
to another — e.g. from [Eq. (6)] to [Eq. (7)]. 

Eor certain choices of exchange parameters, however, 
takes on commensurate values. One particular case, oc¬ 
curring for 


Ji+ 2J2-2J3=0 , [0<2J3 < |Ji| ] , (10) 

[cf. the blue dashed line in Eig. 4], is 


Q 3 s°b= ( ^,0) or (0,^ 
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In this case the ground state of the frustrated ferromag¬ 
net [Eq. (1)] is composed of stripes of spins, with 

three-sublattice order. 

Three-sublattice order also occurs in one of the 
paradigmatic examples of frustrated antiferromagnetism, 
the Heisenberg antiferromagnet on the triangular lattice 

= jY.Si-Sj-hY,st , J>o. (12) 

(u> * 

Eirst studied as a potential route to a quantum spin liq- 
uid^^, this model has a long and distinguished history as 
a testing ground for ideas about classical and quantum 
magnets^’^^. Its behaviour in magnetic field, in partic¬ 
ular, where a collinear 1/3-magnetisation plateau is se¬ 
lected by fluctuations from a degenerate set of classical 
ground states^’^^’^^, has come to be seen as one of the 
paradigms for frustrated magnets. 

At a mean-field level, for a ID spiral state with wave 
vector QsPub [^4- ( 11 ) 1 5 [^4- (1)1 equivalent to 

[Eq. (12)], with an effective interaction 


J = Ja 


Jl + 2J2 + J3 

3 


(13) 


As a consequence, the frustrated ferromagnet, ^ in¬ 
herits many of the special properties of the Heisenberg 
antiferromagnet on a triangular lattice. In particular, the 
classical ground states of 1-L^ can be found by rewriting 
the model as 


A 



h , 


2 


3 


81Ji 


( 11 ) 


(14) 
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where Bz is a unit vector in the direction of the magnetic 
field, the sum runs over all triangles in the lattice, 
and 

niA = ^ X] 

ieA 

is the magnetisation per spin in a given triangle. 

It follows from Eq. (14) that any state for which 

iiiA = ^ Bz V A G lattice , (16) 

9Ja 

is a classical ground state of , and that the classical 
ground state interpolates to saturation (full magnetisa¬ 
tion) for a magnetic field 

hs^t = 9Ja . (17) 

At a mean-field level, both of these properties carry over 
to three-sublattice ground states of [Eq. (1)]. 

However this is by no means the end of the story 
— except at saturation, the mean-field constraint on 
niA [Eq. (16)], does not uniquely constrain the ground 
state, and fluctuations play a crucial role in establish¬ 
ing order. In the case of the Heisenberg antiferromagnet 
on a triangular lattice, [Eq. (12)], it is known that 
thermal^’^^’^^ and/or quantum^^ fluctuations act on the 
large family of degenerate three-sublattice ground states 
to select : 

• a coplanar “ 120 -degree” state with zero magnetisa¬ 
tion, for h = 0. 

• a coplanar “Y-state” with finite magnetisation, in¬ 
terpolating to 120 -degree order for h ^ 0. 

• a collinear “uud” state associated with a 

1/3-magnetisation plateau for intermediate 

values of h. 

• a coplanar 2:1 canted state, interpolating between 
the 1/3-magnetisation plateau and the saturated 
state for h hsat- 

The associated classical phase diagram is shown in Eig. 1 
of Ref. 57. Since strict long-range order of phases which 
break a continuous symmetry is forbidden at finite tem¬ 
perature in two dimensions by the Mermin-Wagner the¬ 
orem, the coplanar phases should be understood as alge¬ 
braically correlated. 

In Sec. HI we use classical Monte Carlo simulation 
to explore how the equivalent three-sublattice order in 
the square-lattice frustrated ferromagnet [Eq. (1)] 

evolves as a function of temperature and magnetic field. 
Eor these purposes, we select the parameter set 

Ja = (Ji, A, Js) = (-1, 3/4, 1/4), (18) 

marked with a black dot in Eig. 4, as representa¬ 
tive of parameters with a three-sublattice ground state 
[cf. Eq. (10)]. This in turn sets a characteristic scale for 
temperature and magnetic field, through Eq. (13). 


C. Highly—degenerate manifold of classical ground 
states 

An enlarged ground state manifold occurs where the 
ID and 2D spirals meet, for 

J2-2J3 = 0 V |Ji| <4J2 <4|Ji| , (19) 

[cf. Eig. 4]. On this phase boundary, the classical ground 
states of 1-L^^ are 2D spirals with wave vector 

= (gr, g™"), 

g™® = ± arccos (^ - cos g™®) . (20) 

which interpolates from the ID to the 2D spiral. This set 
of Q forms a ring in reciprocal space, centred on 

r = (0, 0), (21) 

and defines a highly degenerate manifold of states. 

In Sec. IV we use classical Monte Carlo simulation to 
explore the consequences of this enlarged ground-state 
degeneracy at finite temperature and magnetic field. We 
consider the limiting case where the line Eq. (10), cor¬ 
responding to a ID spiral with three-sublattice order, 
terminates at the phase boundary between ID and 2D 
spirals at 

Jb = (Ji, J 2 , h) = (- 1 , 1 , 1/2) . (22) 

This parameter set is denoted as a black dot in Eig. 4. 
The set of ground-state wave vectors [Eq. (20)] at 
Jb is shown as a dashed line in Eig. 1(b). 

D. Monte Carlo simulation method 

In Sec. HI and Sec. IV of this paper, we use large-scale 
Monte Carlo simulation to study the finite-temperature 
properties of the square lattice frustrated ferromagnet 
'^FFM applied magnetic field. Simulations 

were performed using parallel tempering^^, using 48 to 
80 replicas (temperatures). Simulations were carried out 
for square clusters of 

N = Lx L (23) 

spins, with periodic boundary conditions. The linear 
sizes L were chosen to be commensurate with possible 
Q wave vectors, in the range 60 < L < 180. Typical 
simulations involved 2x10^ steps, half of which were dis¬ 
carded for thermalisation. At every 10 steps there was 
an attempt at exchanging replicas at neighbouring tem¬ 
peratures. Energy scales (field and temperature) were 
normalised to Ja [Eq. (13)], for easy comparison with 
the triangular lattice antiferromagnet^^. 

The simulations for the parameter setJ^ [Eq. (22)], 
presented in Sec. IV, become very challenging at low 
temperatures, especially for low values of magnetic field, 
where several different phases compete. In this case we 
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performed different simulations runs, with and without 
parallel tempering, where the initial state was either one 
of the candidate phases at T = 0, or a configuration com¬ 
posed of domains of different phases, in order to ascertain 
their relative stability. 


E. A short catalogue of order parameters and 
correlation functions 

Experience of simulating the Heisenberg antiferromag- 
net on a triangular lattice [Ref. 57], together with the 
symmetry of the ordered phases found in the absence of 
magnetic field [cf. Fig. 4], suggests a number of order 
parameters and correlation functions likely to be of use 
in determining the behaviour of 1-L^^ [Eq. (1)] in applied 
magnetic field. 

The ID spiral phase in Fig. 4 breaks the four-fold ro¬ 
tational symmetry of the square lattice from C 4 down to 
C 2 , by choosing to orient the stripes in vertical or hori¬ 
zontal direction. We analyse this with the following order 
parameter 


Oc. = (|<^x|) , (24) 

'^x=^ESi-(Si+x-Si+y). (25) 

i 

The sum over % runs over all N = lattice sites and 
the lattice unit vectors are denoted as x = ( 1 , 0 ) and 

y = (0,1). 

In the presence of field, the three-sublattice ordered 
states break the translational symmetry of the lattice 
along the direction. This can be measured by an order 
parameter based on a two-dimensional irreducible repre¬ 
sentation of (73 = ^ 3 , cf. the triangular-lattice case^^ 


^x,l 


(26) 

II 

CM 


(27) 


where the sum over i runs over the N/3 “3-spin stripes”, 
each equivalent to an elementary triangle. In order to 
account for the two ways of breaking the C 2 symmetry, 
the Z 3 order parameter must be measured along both 
directions of the lattice, i.e. x —y in Eqs. (26) and (27). 
The resulting final order parameter is 




KlP + l^: 


x,2l 


IV’y.ll 


iv-. 


-, 1 / 2 ' 


y,2l 


(28) 


The susceptibility associated with each order parameter 
O is defined as 


X = N 


{m-{\o\Y 


(29) 


In the presence of magnetic field, quasi-long-range or¬ 
der can develop in the transverse components of spin 

si = {8^,81). (30) 


This is captured by the spin stiffness see e.g. Refs. 57 
and 62 and references therein. 


Ps[e\ 


= v(E ^ (g ■ S 




\ <5 (ij), / 

which is averaged over the two lattice directions 

e = {x,y}. 

The specific heat, defined as 


Ch = 


1 (E) - {Ej 

N r2 


(32) 


where E is the total internal energy, is also helpful in 
tracking phase transitions. 

Last, but not least, the momentum-resolved spin struc¬ 
ture factor 


^(q) = (;^ I E exp(-zq • 

i 



(33) 


is of considerable importance for characterising the spin- 
liquid phase studied in Sec. IV. Where appropriate, we 
analyse separately the structure factors for the longitu¬ 
dinal and transverse components of spin 


*5^ (q) = ( ^ I E exp(-iq ■ Tj) 

i 

' 5 -^(q) = (;^|E®texp(-iq-ri) 


(34) 

(35) 


where is defined in Eq. (30). 


III. TRIANGULAR-LATTICE PHYSICS ON A 
SQUARE LATTICE 

In Sec. IIB, we established a connection between the 
square-lattice frustrated ferromagnet, [Eq. (I)], 

and the Heisenberg antiferromagnet on a triangular lat¬ 
tice, [Eq. (12)1, in the case where the ground state 
of is a ID-spiral with wave vector 

QSb=(f.o) or , 

[cf. Eq. (II)], corresponding to three-sublattice “stripe” 
order. An example of a three-sublattice stripe state, with 
finite magnetisation, is shown in Fig. 1(a). The momen¬ 
tum set of the ground state manifold 

Qa = {(±2^/3,0), (0,±2^/3)}, 


T 


(36) 
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FIG. 5. (Color online). Finite-size phase diagram of a 
square-lattice frustrated ferromagnet in applied magnetic 
field. The phases found — a coplanar “Y-state”; a collinear 
1/3-magnetisation plateau; and a coplanar 2:1 canted state 
— and the structure of the phase diagram, closely parallel 
finite-size results for the Heisenberg antiferromagnet on a tri¬ 
angular lattice^^. In both cases, finite-size effects strongly 
renormalise the temperature associated with the transition 
from the collinear 1/3-magnetisation plateau into the Y- 
state. Results are taken from classical Monte Carlo simu¬ 
lations of [Eq. (1)] for the parameter set Ja [Eq. (18)]. 

Phase boundaries were extracted from anomalies in the or¬ 
der parameter associated with lattice rotations [Eq. (24)] and 
spin stiffness ps [Eq. (30)], for a cluster of Y = 90^ = 8100 
spins. Temperature and magnetic field are measured in units 
of Ja [Eq. (13)]. The horizontal dashed line corresponds to 
the temperature cut used in Eig. 7. 


has four components, whereas that of the triangular an¬ 
tiferromagnet has only two. 

In what follows, we use classical Monte Carlo simu¬ 
lation to explore the properties of at finite tem¬ 

perature and magnetic field, considering a parameter 
set Ja [Eq. (18)]. The results of these simulations are 
summarised in the phase diagram Fig. 2. The similarities 
to the magnetic phase diagram of the triangular-lattice 
antiferromagnet cf. Fig. 1 of Ref. 57, are striking. 
At first sight, the main difference is only that the order¬ 
ing temperature scale is roughly double that for . 

As in the case of the triangular-lattice antiferromagnet^^, 
it will be instructive to compare the phase diagram where 



FIG. 6. (Color online). Finite-temperature phase tran¬ 
sition into a three-sublattice, “120-degree” ground state, 
breaking lattice rotation symmetry. The order parameter 
Oc 2 [Eq. (24)] shows finite-size scaling consistent with a 
phase transition in the Ising universality class, in contrast 
with the triangular-lattice Heisenberg modeE^. Results are 
taken from classical Monte Carlo simulations of [Eq. (1)] 
for the parameter set Ja [Eq. (18)], in the absence of magnetic 
field (/i = 0). 


a proper finite-size scaling has been performed [Fig. 2], 
to one extracted from simulations of a fixed cluster size 
L = 90 [Fig. 5]. 

While the two models have much in common, differ¬ 
ences arise in the way in which ordered phases break 
lattice symmetries. In both cases, under applied mag¬ 
netic field, three-sublattice ordered phases break a dis¬ 
crete Cs = Z 3 symmetry, associated with the interchange 
of the different sublattices. However the stripe-like or¬ 
der found in the square-lattice model also breaks a C 2 
lattice-rotation symmetry, when choosing between the 
two possible ordering vectors, [Eq* (H)]- This 

additional symmetry has a number of interesting con¬ 
sequences, described below. 


A. Ising transition at /i = 0 

We consider first the limit of vanishing magnetic field, 
corresponding to /i/Ja = 0 in Fig. 2. In the absence of 
magnetic field, the Heisenberg antiferromagnet on a tri¬ 
angular lattice [Eq. (12)] has been argued to exhibit 
a finite-temperature transition, linked with the prolif¬ 
eration of Z 2 vortices associated with spin chirality^^. 
Whether this process corresponds to either a true phase 
transition or a crossover is still under debate^^“^^. 

The situation in the square-lattice frustrated ferro¬ 
magnet [Eq. ( 1 )] is quite different. At = 0 , a clear 

phase transition is observed, associated with the breaking 
of lattice-rotation symmetry by three-sublattice stripe 
order. The relevant order parameter is Oc 2 [Eq. (24)], 
and in Fig. 6 we show a scaling plot of simulation results 









T/Ja 


(c) 



(d) E/L^ 



0 0.04 0.08 0.12 
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FIG. 7. (Color online). Evidence for a collinear 1/3-magnetisation plateau and coplanar “Y-state” in a square-lattice frustrated 
ferromagnet. (a) Temperature dependence of the order parameters Oc 2 [Eq. (24)] and [Eq. (28)], showing the onset of 
a collinear state with three-sublattice stripe order (1/3-magnetisation plateau), at T 0.78Ja- (b) Temperature dependence 
of the spin stiffness ps [Eq. (31)], consistent with a BKT transition into the “Y-state” at Tl ~ 0.24Ja. ( c ) Energy histogram 
calculated for T = 0.776Ja, showing the discontinuous nature of the phase transition from paramagnet to 1/3-magnetisation 
plateau, (d) Einite-size scaling of the temperature associated with the BKT transition, Tl [Eq. (41)]. Results are taken from 
classical Monte Carlo simulations of 1-i^^ [Eq. (1)] for the parameter set Ja [Eq. (18)], in applied magnetic field /i/Ja = 2, 
cf. dashed line in Eig. 5. Temperature is measured in units of Ja [Eq. (13)]. 


for 

Oc^ = (37) 

as a function of the reduced temperature 



where the Ising critical exponents 

= 1 , ^ ^ 

are found to describe the data well. 

We do not observe behaviour explicitly related to the 
unbinding of Z 2 vortices in simulation of TL^^. However, 
just as in simulations of TL^ [cf. Ref. 57], the correlation 
length associated with transverse components of spin be¬ 
comes very large as ^ 0 , making it difficult to draw 
definitive conclusions. 


B. Double transition for 0 < /i < 3 Ja 

We now consider the phases found for low to intermedi¬ 
ate values of magnetic field in the phase diagram Fig. 2. 
Here, a double phase transition is observed as the sys¬ 
tem is cooled down from the paramagnet — cf. Fig. 7, 


for 0 < /i/ Ja <3. As the system is cooled from the high- 
temperature paramagnet, it first enters the collinear 1/3- 
magnetisation plateau, as described by the rise of the or¬ 
der parameters in Fig. 7(a). This process is a single phase 
transition, which simultaneously breaks the C 2 lattice- 
rotational symmetry [Eq. (24)], and the Z 3 translational 
symmetry [Eq. (28)] associated with the three-sublattice 
order. For values of magnetic field 1.5 ^ /?^/Ja < 3, a 
double distribution in the internal energy is clearly ob¬ 
served at the transition temperature [Fig. 7(c)], which 
shows that the transition is of first order. We cannot 
observe this behaviour for lower magnetic fields, presum¬ 
ably due to the increased finite-size effects. 

In the case of the Heisenberg antiferromagnet on a 
triangular lattice, different 1/3-magnetisation plateau 
states are connected by a 3-fold permutation symme¬ 
try, and the phase transition into the paramagnet is con¬ 
tinuous^^. Meanwhile, in the case of the square-lattice 
frustrated ferromagnet, the permutation of different sub¬ 
lattices is complemented by a lattice rotation, enlarging 
the symmetry from Z 3 to Z 2 x Z 3 , and the correspond¬ 
ing phase transition is first order. This is reminiscent 
of the 6 -state Potts model in 2D, whose ordering phase 
transition is known to be first order^^’^^. 

The canted Y state is found by lowering the temper¬ 
ature further from the plateau phase. In addition of 
breaking the C 2 and Z 3 symmetries, this phase displays 
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algebraic order in the S^-S^ plane, as shown by the fi¬ 
nite spin stiffness in Fig. 7(b). Assuming a Berezinski- 
Kosterlitz-Thouless (BKT) transition^^, the transition 
temperature in a finite-size system is found via the jump 
in spin stiffness 

Ti = I X Aps, (40) 

which can then be finite-size scaled in a characteristic 
logarithmic fashion^^’^^ 

TL=TBKT{l + l^ogLl\ogb)^ (41) 

as shown in Fig. 7(d). 

The finite-size corrections to this transition are ob¬ 
served to be rather large. From simulations of a fixed 
cluster size L, cf. the phase diagram for L=90 in Fig. 5, it 
is unclear if the BKT transition merges with the C 2 ^3 
transition for fields h < 0.6, i.e. if there is a single 
transition from the paramagnet to the Y state, or two. 
However, these two phase transitions are observed to be 
clearly separated once the proper finite-size scalings are 
performed, and no direct transition from the paramagnet 
into the Y phase can be reported for fields h > 0.4Ja 
[Fig. 2]. We have faded the region in Fig. 2 correspond¬ 
ing to values of field 0 < h < 0.4Ja, since finite-size 
effects become very large in this region and this analysis 
becomes unreliable. However, the available data still sug¬ 
gests the presence of a double phase transition as the field 
approaches h ^ 0. This behaviour very closely matches 
what is observed in the Heisenberg antiferromagnet on 
a triangular lattice, which has been discussed in detail 
previously^^. 

For magnetic fields above the plateau h > 3Ja, we reg¬ 
ister a single phase transition into the 2:1 canted state. 
The location of several quantities, such as the jump in 
the spin stiffness, anomalies in the C 2 and Z3 order- 
parameter susceptibilities, and in the specific heat, coin¬ 
cide within resolution, after they are properly finite-size 
scaled. We do not observe any discontinuity in the in¬ 
ternal energy as a function of temperature, or a bimodal 
energy distribution at the transition temperature, in the 
clusters simulated. None the less, the maximum value 
of the specific heat scales roughly with and the tem¬ 
perature at which it is found scales roughly with 1/A^, 
behaviour consistent with a first-order transition. 


IV. CONSEQUENCES OF AN ENLARGED 
GROUND STATE MANIFOLD 

We now consider the properties of 1-L^^ [Eq. (1)] for 
the parameter set [Eq. (22)], as summarised in the 
phase diagram, Eig. 8. This parameter set corresponds 
to the point in the classical ground state phase diagram 
[Eig. 4], where the line of ID spirals with three-sublattice 
order terminates on the boundary with the 2D spiral 
phase. The behaviour of the model at this special point is 



FIG. 8. (Color online). Finite-size phase diagram of 
a square-lattice frustrated ferromagnet in applied magnetic 
field, showing a classical spin liquid and, at low fields, a vortex 
crystal [cf. Fig. 1(c)]. The other phases found are a conical spi¬ 
ral state interpolating to zero magnetisation; a coplanar “Y- 
state” and a collinear 1/3-magnetisation plateau [cf. Fig. 1(a)] 
at intermediate magnetisation; and a coplanar 2:1 canted 
state interpolating to saturation. Results are taken from clas¬ 
sical Monte Carlo simulations of 1-L^^ [Eq. (1)] for the param¬ 
eter set J B [Eq. (22)], corresponding to a point of high degen¬ 
eracy. Phase boundaries were extracted from anomalies in the 
order parameter associated with lattice rotations [Eq. (24)], 
spin stiffness ps [Eq. (30)] and the spin structure factor 
<S(q) [Eq. (33)], for a cluster of iV = = 120^ = 14400 

spins. Temperature and magnetic field are measured in units 
of Ja [Eq. (13)]. The horizontal dashed lines corresponds to 
the temperature cuts used in Figs. 9-13 and Fig. 14. 

determined by the highly-degenerate manifold of states 
described in Sec. HC. And, as we shall see, the “ring” 
structure defined in Eq. (20) leaves a characteristic fin¬ 
gerprint in the spin structure factor, for all temperatures 
and values of magnetic field. 

Eor a generic choice of parameters on the boundary 
between ID and 2D spirals, it is not possible to find 
commensurate wave vectors Q which satisfy Eq. (20). 
However for specific choices of parameters J, commensu¬ 
rate wave vectors do exist, and the parameter set which 
offers the greatest number of commensurate states is in 
fact set Jb- Here the wave vectors 

Qr = {(0,27r/3), (7r/2,7r/3), (tt/S, 37r/5)}, (42) 
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together with their reflections and mirrors, fulfill 
Eq. (20). In total, this give 20 different commensurate 
wave vectors, compatible with a wide range of cluster 
sizes. 

We begin by establishing the finite-temperature phase 
diagram. Fig. 8, using Monte Carlo simulations for a clus¬ 
ter of size L = 120. For values of applied field h > 1.2Ja, 
the phase diagram closely follows that of the parameter 
set Ja [Eq. (18)] (and therefore the Heisenberg antifer- 
romagnet on a triangular lattice, [Eq. (12)]), dis¬ 

playing the succession of canted Y phase, collinear 1/3- 
magnetisation plateau and canted 2 : 1 phase as a func¬ 
tion of increasing magnetic field. 

Simulations become very challenging at lower values 
of applied field, h < 1.2Ja, making it very difficult to 
extrapolate results to the thermodynamic limit. Still, in 
this region we can establish the presence of two additional 
non-coplanar phases. At very low temperatures, a conical 
spiral phase described by a single wave-vector is favoured 
by thermal fluctuations, as described in Section IV B. At 
intermediate temperatures, between this phase and the 
disordered paramagnet phase, a novel conical phase, with 
a more complex spin texture, is stabilised, as described 
in Section IV C — cf. Fig. 1(c). 

All the ordered phases at are described by wave 
vectors belonging to, or combinations of, [Eq. (42)]. 
All the transitions from the disordered phase are very 
clearly observed to be of first order. The overall ordering- 
temperature scale is strongly reduced from the case, 
both in units of Ja and |Ji|, indicating the larger role 
played by frustration in this case, suppressing the order¬ 
ing tendency of the model. 


A. Spin-liquid phase from a ring of correlations 

We consider first the nature of the paramagnetic phase 
found for all values of magnetic field at sufficiently high 
temperature, as shown in the phase diagram Fig. 8. 
The highly-degenerate ground-state manifold of 1-L^^ at 
the parameter set [Eq. (22)] is not just of concern at 
T = 0, but also has profound consequences at finite tem¬ 
perature. Its effects are most obvious in the spin struc¬ 
ture factor S{q) [Eq. (33)] which takes on a finite value 
for all wave vectors q which satisfy, or nearly satisfy, the 
“ring” condition Eq. (22). A corresponding ring structure 
can be seen in S{q) for all of the ordered phases shown 
in Fig. 8, and is equally prominent for h = 0 [Fig. 1(b)], 
and for h = 2.4Ja [Fig. 9 . In the paramagnetic regions 
of the phase diagram. Fig. 8, where no symmetries are 
broken, the extra degeneracy gives rise to a classical spin 
liquid. 

The nature of the correlations along the ring does not 
change in any fundamental way when magnetic field is 
applied. In Fig. 9 we show the evolution of the full 
structure factor S{q) at h = 2.4Ja, as the tempera¬ 
ture is lowered from the disordered phase into the 1/3- 
magnetisation plateau. At high temperatures and away 


from the transition [Fig. 9(a)], the height of 5(q) is very 
uniform, such that many wave vectors belonging to the 
ring, or very close to it, contribute equally. If the tem¬ 
perature is raised from this point, the height of S{q) 
just decreases smoothly, and no crossover into a standard 
paramagnet is observed at any temperature scale. When 
the system is cooled down to near the phase transition 
T ^ 0.29Ja, Fig. 9(b), peaks begin to develop at the 
wave vectors Qa [Eq. (36)], corresponding to incipient 
three-sublattice magnetic order. Inside the magnetisa¬ 
tion plateau at even lower temperature. Fig. 9(c), Bragg 
peaks (from the (q) component of the structure factor) 
are observed, associated with the broken Z3 symmetry 
discussed in Sec. III. A diffuse ring of low-lying excita¬ 
tions can still be observed in this phase, attesting to the 
pervasiveness of the ring correlations. 

In order to further clarify the physical implications of 
this ring in the structure factor, we now focus on the 
transverse component iS^(q) [Eq. (35)]. Fig. 10(a) shows 
a cut of 5^(q) along the line = 0 in the disordered 
phase for h = 2.4Ja and T = 0.45Ja- We find that 
the structure factor can be fitted quite well to a (double) 
Lorentzian expression of the form 




+ {q - qoY' 


(43) 


centred at qo = ±27r/3. Similar fits can be performed for 
different cuts of iS^(q) in the Qx — Qy plane. 

The width of the structure factor around go is con¬ 
trolled by the correlation length which sets a finite 
length-scale for correlations. S^{q) shows practically no 
finite-size dependence for the clusters studied, attesting 
to the short-range nature of the correlations. The correla¬ 
tion length ^ decays slowly with temperature [Fig. 10(b)], 
which indicates the stability of the ring correlations for 
a wide range of temperatures. 

Our next target is the total structure factor arising 
from low-energy excitations associated with the ring. 
First, we need identify the wave vectors Qring contribut¬ 
ing to the total ring structure factor. A finite ^ means 
that there are momenta outside the T=0 ring [Eq. (20)] 
which will contribute to the total structure factor, [cf. 
Fig. 1(b)]. In order to capture these wave vectors with 
a finite spectral weight, we take the one-sublattice spin- 
wave dispersion Eq. (A57) [see Appendix A] and impose 
an energy cutoff of ujc = O.OIJa on it. Alternatively, 
using the inverse of the correlation length of Fig. 10(b) 
as input, we can predict approximately the wave-vectors 
which will be thermally activated. Both approaches 
provide the same set of wave vectors belonging to the 
enlarged ring, but the latter, correlation-length based 
method becomes unwieldy and computationally more ex¬ 
pensive for large system sizes. 

The number of momentum vectors, thermally ac¬ 
tivated is found to scale linearly with the cluster size 
Nr ^ for a very wide range of system sizes [Fig. 11]. 
The area covered by the ring is given roughly by the prod¬ 
uct of its width with its perimeter, which scales with L. 
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FIG. 9. (Color online). Transition from the helicoidal spin liquid into a collinear m = 1/3 magnetisation plateau, in applied 
magnetic field, as revealed by the spin structure factor <S(q) [Eq. (33)]. (a) Spin correlations for T = 0.44Ja, showing a “ring” 
feature characteristic of the high-temperature spin liquid, (b) Spin correlations for T = 0.29Ja, approaching the transition 
into an ordered state, showing an enhancement of fluctuations near to potential ordering wave vectors, (c) Spin correlations 
for T = 0.24Ja, within the collinear m = 1/3 magnetisation plateau, showing Bragg peaks at Qa^ub Pq- (H)], coexisting with 
the “ring” feature. Results are taken from classical Monte Carlo simulations of 1-L^^ [Eq. (1)], for a cluster of — 90^ spins, 
with exchange parameters Js [Eq. (22)]. Magnetic field was set to h — 2.4Ja, corresponding to the dashed in line in Eig. 8. 
The q = 0 component of the spin correlations has been subtracted for clarity. 


(a) (b) 




EIG. 10. (Color online) Characterisation of the helicoidal spin liquid in applied magnetic field h = 2.4Ja. (a) Structure factor 
‘^^(^x,0) [Eq. (35)], showing a cross section of the characteristic “ring” of scattering [cf. Eig. 9(a)]. Data for a wide range of 
system sizes collapses onto a single curve. Dashed lines are fits to a Lorentzian [Eq. (43)], from which the correlation length ^ 
is obtained, (b) Temperature-dependence of the correlation length Results are taken from classical Monte Carlo simulations 
of [Eq. (1)], for clusters of iV = L x L spins, with T = 0.45Ja and the same values of exchange and magnetic field as 

Eig. 9. 


This implies that there are approximately L cuts across 
the ring such as Fig. 10(a). The ring has a fixed width in 
momentum space, set by which is basically indepen¬ 
dent of L, see Fig. 10(b). Since the resolution in momen¬ 
tum space scales with system size as 1/“^, the number of 
points allowed inside each cut of the form Fig. 10(a) must 
scale as L. Therefore, the number of points covered by 
the ring scales as Nr r-^ L x L, and is an extensive prop¬ 
erty of the system. 

In order to gain more insight into this state, we de¬ 
fine an integral of the transverse structure factor 5'^(q) 


[Eq. (35)] over the manifold of states contributing to the 
ring 

S^ = ^<S-^(q) V q e {q|a;isL(q) < Wc} , ( 44 ) 

q 

as well as a sum over the discrete wave vectors Qa 
[E q. (36)] associated with three-sublattice order 

= El 

qCfQA} 
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and the difference of the two 

SLa = - si . (46) 

Fig. 12 shows how these measures of the ring, normalised 
to the number of contributing wave vectors, evolve across 
the transition from the high temperature spin-liquid into 
the 1/3-magnetisation plateau for h = 2.4Ja. It is clear 
that the 5(q) ring has a uniform height at temperatures 
away from the phase transition. Right above the tran¬ 
sition, enhanced fluctuations in herald three- 

sublattice order, but the defining characteristics of the 
ring survive both in the disordered and ordered limits. 

The structure factor per spin 

/Nr ~ Sfl/A (47) 

evaluated at a typical wdive vector Qring [Eq. (20)], is itself 
not an extensive quantity. Instead it vanishes as l/T^, 
as can be seen in Fig. 13(a) for a variety of temperatures 
in the disordered region. This is also the behaviour ob¬ 
served in a standard paramagnetic region. However, the 
total ring structure factor per spin is a function of 

Nji ^ and should therefore be an extensive quantity. 
Fig. 13(b) shows that converges to a finite value 

as 1/“^ ^ 0, and therefore 'Er is a non-zero quantity in 
the thermodynamic limit. This is a direct demonstration 
of an extensive, non-trivial degeneracy in the thermody¬ 
namic limit, an hallmark of a spin liquid. 


B. Low—temperature conical spiral 

We now turn to the nature of the low-temperature 
ordered state found for h/J^ < 1 and T/J^ < 1.7 in 
the phase diagram Fig. 8. We find that specific ordered 



FIG. 11. (Color online) Evidence that the number of spin 
configurations contributing to the helicoidal spin liquid, Nr, 
scales linearly with the system-size, N — . Nr was esti¬ 

mated on the basis of the number of states lying close to the 
ring defined by [Eq. (20)], for parameters Js [Eq. (22)] and 
system sizes ranging from L = 12^ to L = 450^, as described 
in Section IV A. Inset : wave vectors of states contributing 
to Nr for L = 60. 



EIG. 12. (Color online). Comparison between the spec¬ 
tral weight associated with the ordering vectors of the 1/3- 
magnetisation plateau, and that associated with the remain¬ 
ing wave-vectors belonging to the ring [cf. Eig. 9]. The total 
spectral weight in the ring, [Eq. (44)], undergoes a dis¬ 
continuous drop at the ordering temperature T = 2.8Ja, but 
remains finite in the ordered state. Eluctuations at the or¬ 
dering vector, Er_^ [Eq. (46)], are associated with a discon¬ 
tinuous rise at the ordering temperature. Results are taken 
from classical Monte Carlo simulations of [Eq. (1)] for 

system sizes L = 48, 60, 70, and parameters [Eq. (22)] 
with h = 2.4Ja. 

states, with a wave vector belonging to the ring-manifold, 
Eq. 20, are selected at low temperature. None the less, 
“ring” seen in the spin structure factor S{q) [Eq. (33)] 
survives in these ordered phases, as shown in Fig. 12. For 
magnetic fields h > 1.2Ja, the magnetisation process of 
the triangular lattice antiferromagnet is recovered, as de¬ 
scribed in Section III for the parameter set [Eq. (18)]. 
For values of field 0 < h < 1.2Ja, we observe a strong 
competition between the coplanar Y state and different 
conical states, all with uniform magnetisation Sf. 

In order to better understand this behaviour, we have 
carried out a low-temperature expansion of the spin-wave 
excitations around the T = 0 classical ground-state man¬ 
ifold, which allows the calculation of corrections to the 
free energy as a power series in T. Details of these calcu¬ 
lations are given in Appendix A 1. In the absence of mag¬ 
netic field, this expansion predicts that low-temperature 
fluctuations favour a coplanar spiral state described by a 
single incommensurate wave vector very close to 

QiQ = (7r/3,7r/2) (48) 

(plus associated mirrors and reflections), cf. Fig 15. 

When a small magnetic field is applied, the state 
favoured by fluctuations in the harmonic approach is 
a conical, non-coplanar version of the h = 0 state, cf. 
Fig. 14(a) and (d). The wave-vector of the h = 0 spiral 
is preserved in the plane, while all spins have a 

constant magnetisation along the spin direction. This 
state has been identified previously as the conical um¬ 
brella state in the anisotropic triangular lattice^’^^’^'^’^^. 
According to the low-temperature expansion, the entropy 






13 



FIG. 13. (Color online). Evidence for spin-liquid behaviour 
in the finite-size scaling of the integrated structure factor, 
[Eq. (44)]. (a) The total weight associated with spin 

configurations contributing to the “ring” vanishes as 1/L^ in 
the thermodynamic limit, (b) The structure factor per spin 
summed over all ring points, E;^/L^, scales linearly with the 
system size and thus persists in the thermodynamic limit. 
Lines are a linear fit to the data for the largest system sizes. 
Results are taken from classical Monte Carlo simulations of 
1-L^^ [Eq. (1)1 for h = 2.4Ja and T = 0.45Ja, with exchange 
parameters [Eq. (22)]. 


of the Y coplanar state increases with increasing value of 
magnetic field, while that of the conical state decreases. 
For fields h > O.SJa, the Y state is predicted to be the 
preferred low-temperature state. 

Monte Carlo simulations confirm that, for values of 
field 0 < < O.SJa, the conical state is stabilised at 

low temperature, while for 0.8 ^ h < 3 Ja the Y state is 
favoured. The first-order phase transition between both 
phases is very difficult to observe in MC simulations, 
since it happens after two other symmetry-breaking 
phase transitions at higher temperature. We also find, 
from e.g. a negative spin stiffness (not shown), that simu¬ 
lations seeded with a conical spin texture with an in-plane 
wave vector Qig = (7r/3,7r/2) show a strong tendency to 
become slightly incommensurate in the S^—S^ plane, in 
full agreement with the harmonic approximation. This 
selection of an incommensurate, non-coplanar phase by 
thermal fluctuations is an interesting counterexample to 
the rule of thumb that fluctuations prefer collinear, or 
at worst coplanar, phases, as a manifestation of the cel¬ 
ebrated order-due-to-disorder effect, cf. Ref. 59. 


C. Magnetic vortex crystal 

We next consider the nature of the ordered state found 
for intermediate temperatures, 1.7 < T/Ja ^ 2.2, and 


low values of magnetic field, h/J/\ < 1.2, between the 
low-temperature conical spiral and the high-temperature 
spin-liquid phase, as shown in the phase diagram Fig. 8. 
In this parameter range, Monte Carlo simulations started 
from a random initial configuration often get trapped 
in local free-energy minima, resulting in domain walls 
between different phases. This hints at the presence of 
competing phases, and a first-order transition, at a lower 
temperature than the initial ordering transition. Fur¬ 
ther simulations with open boundary conditions'^ reveal 
a tendency for the edge spins to be collinear with 5'^, 
while the bulk reproduces the same behaviour as for peri¬ 
odic boundary conditions. We resort to comparing differ¬ 
ent parallel-tempering simulations for each value of field, 
initialised from a variety of different ordered states, in¬ 
cluding random configurations, mixing different states in 
a single replica, or across temperature space. This ap¬ 
proach does not permit a very accurate location of the 
phase transition, but gives us confidence on the phases 
present and on the overall topology of the phase diagram. 

Carrying out this analysis, we find a constant- 
magnetisation state with a spin texture in the S^-S'^ 
plane described by fourwsive vectors, cf. Fig. 1(c). Usu¬ 
ally, multiple-Q states violate the fixed spin-length con¬ 
straint jS^I = 1, and therefore are not favoured at low 
temperatures. However, it is possible to construct a very 
specific 4(5 state out of spirals with the wave vector 

Q 4 Q = (37r/5,7r/5) (49) 

and symmetry-related vectors. For a given magnetisa¬ 
tion 


m = Sf = const. , (50) 

the spin configuration associated with this AQ state is 

^2^"Re[F(ri)] , (51) 

^"^m[F(r,)], (52) 


S! = 
Sf = 


1 — 

l-rri^ 


where 


F(rj)= + e*2^/5g-i.QA.ri 

_^gi47r/15gi.QB.ri _^g-i87r/15g-iQB.ri ^ ( 53 ) 


and 


Qa = (^/5, -3^/5) , 

Qb = (3^/5,V5). (54) 

For vanishing magnetisation {m = Sf = 0), the AQ state 
has the form of a lattice of vortices shown in Fig. 14(e). 
By construction, the AQ state interpolates smoothly to fi¬ 
nite magnetisation, where it takes on the form illustrated 
in Fig. 1(c). The associated spin structure factor has four 
peaks [cf. Fig. 14(b)]. 

Monte Carlo simulations for a range of fields 
0 < < 1.2Ja, initialised with half of the replicas in the 

IQ state, and the other half in the AQ state, paint a con¬ 
sistent picture of a double phase transition, as shown in 
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FIG. 14. (Color online). Finite-temperature phase transition between conical spiral state and magnetic vortex crystal found 
for low value of magnetic field in the phase diagram Fig. 8. (a) Structure factor <S^(q) [Eq. (35)], within the conical-spiral state 
for T = 0.13 Ja, showing peaks at Qiq [Eq. (48)]. (b) Structure factor within the magnetic vortex crystal for T = 0.175 Ja, 
showing Bragg peaks at Q 4 Q [Eq. (49)]. (c) Temperature-dependence of the peaks, showing a phase transition from conical 
spiral to magnetic vortex crystal at T 1.6 Ja [of. Eig 8]. (d) Illustrative spin configuration in the plane perpendicular 
to magnetic field associated with the IQ conical-spiral state, (e) Illustrative spin configuration in the plane perpendicular 
to magnetic field associated with the 4Q magnetic vortex crystal. Results for <S"'"(q) are taken from classical Monte Carlo 
simulations of 1-i^^ [Eq. (1)] for h = 2.4Ja, with exchange parameters Jb [Eq. (22)]. 


Fig. 14(c), for h = O.SJa- As the system is cooled down 
from the paramagnet, a region is found where the AQ 
state is stabilised, followed by the IQ state at a lower 
temperature. We conclude that the inner phase transi¬ 
tion is probably weakly first-order, since we observe no 
discontinuities in the internal energy, nor any significant 
release of entropy. Unlike simulations of the IQ state 
initialised with the Qig = (7r/3,7r/2) wave vector, the 
AQ state does not show a tendency to become incom¬ 
mensurate, which attests to its stability. We observe this 
scenario of two ordered phases down to the = 0 limit, 
where both phases lose their conical character by becom¬ 
ing coplanar. 

The spin textures of the two non-coplanar states are 
fundamentally different. A spiral phase described by a 
single wave vector breaks a reflection symmetry, which 
can be detected, for example, by defining a chiral order 
parameter on the x-bonds 

Ktotal = ^ E 

i 

Conversely, multiple-Q states typically have vanishing 
net chirality (/^totai) = 0- The AQ state is no exception, 
and the chirality defined on each plaquette points along 


the direction but with alternating signs, so that the 
total chirality is zero. This describes a spin texture where 
the spin texture winds up in alternating directions across 
the lattice, which may be visualised as a crystalline struc¬ 
ture of alternating magnetic vortices [Fig. 14(e)]. We 
have not found other possible multiple-Q states with uni¬ 
form which could be stabilised in this model. 


V. DISCUSSION OF RESULTS 

In this Article, we have explored some of the novel 
phases which can arise in a square-lattice frustrated fer- 
romagnet in applied magnetic field. In the process, we 
have seen a triangular—lattice antiferromagnet reborn on 
a square lattice, and uncovered a spin-liquid and a vortex 
crystal, both with finite magnetisation. In what follows, 
we examine how these phases compare with known re¬ 
sults for frustrated antiferromagnets, and sketch some of 
the possible consequences of quantum fluctuations, pre¬ 
viously touched upon in Refs. 45 and 46. 
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A. Three-sublattice physics and 1/3—magnetisation 
plateau 

The phase diagram for the square-lattice frustrated 
ferromagnet 1-L^^ [Eq. (1)], for the parameter set 
[Eq. (18)], shown in Eig. 2, is markedly different from 
what might normally be expected in a square-lattice frus¬ 
trated magnetHowever, it is remarkably similar to 
that of the triangular-lattice antiferromagnet, developed 
in Ref. 57. This similarity extends beyond the phases 
found — canted “Y” state, collinear 1/3-magnetisation 
plateau, and 2 : 1 canted phase — and the topology of 
the resulting phase diagram, into the finite-size scaling 
of the associated phase transitions (compare Eig. 2 and 
Eig. 5 of this Article with Eig. 1 and Eig. 2 of Ref. 57). 

These results are striking, since the collinear 1/3- 
magnetisation plateau is usually thought to be a hall¬ 
mark of triangular-lattice physics, stabilised by a delicate 
order-by-disorder effect^’^^’^^. Here, however, a robust 
1/3-magnetisation plateau, with three-sublattice order, 
arises in a model more naturally associated with incom¬ 
mensurate phases. In fact, for this particular parameter 
set, the only significant difference between the square- 
lattice and triangular lattice models is the way in which 
ordered phases break lattice-rotation symmetry. 

Given the degree of fine-tuning in the choice of 
parameters, it is natural to ask what happens when 
^FFM ^ 2)1 is tuned away from the special line 

J?, = J 2 — 1 -^ 1 1/2 [cf. Eig. 4]. In this case, the collinear 
and coplanar phases of Eig. 2 would no longer belong to 
the T = 0 ground-state manifold. However, commensu¬ 
rate, collinear phases with Q = (^,0) or (0, ^) are still 
generically more favoured by thermal fluctuations than 
non-collinear or non-coplanar phases^^. And in the case 
of the Heisenberg antiferromagnet on an anisotropic tri¬ 
angular lattice, it is known that quantum fluctuations 
can stabilise a 1/3-magnetisation plateau, even where it 
is not a classical ground state^^’^^’^^. 

In the present case, our preliminary Monte Carlo sim¬ 
ulations hint at a rich phase diagram near the line 
-^3 = -^2 — 1 -^ 1 1/2. At the lowest temperatures, coni¬ 
cal phases with incommensurate wave-vector predomi¬ 
nate. However at higher temperatures, commensurate, 
coplanar phases, and in particular the collinear 1/3- 
magnetisation plateau, are restored by fluctuations^^. 
The full determination of the corresponding phase dia¬ 
gram remains an open problem, and may require use of 
an algorithm specially adapted to incommensurability^^. 

Another question of obvious interest is what hap¬ 
pens for quantum spins. Eor large S', it seems rea¬ 
sonable to expect that quantum fluctuations would act 
much like thermal fluctuations, stabilising “triangular- 
lattice” physics in the square-lattice frustrated ferromag¬ 
net, both on the special line J 3 = J 2 — |^i|/ 2 , and near 
to it. The extreme quantum limit of 1-L^^ [Eq. (1)], 
for spin S = 1 / 2 , has already been studied in exact di- 
agonalisation"^^. These studies confirm the existence of 
1/3-magnetisation plateau at T = 0, for intermediate 


values of magnetic field, and a range of parameters close 
to the line J 3 = J 2 — |Ji|/2. However, at higher values 
of field, instead of interpolating to saturation through a 
canted phase, the magnetisation plateau gives way to a 
quantum spin-nematic phase [cf. Ref. 9 and 46]. 

It is also reasonable to ask what phases arise for 
other parameter sets where the classical ground state of 
^^FM commensurate with the lattice. These 

will occur on lines 

Ji -\- 2J2 + 4 J 3 cos ^^ = 0 , (^ 6 ) 

for all integer p > 3, with associated ordering vector 

Qp=(y,o). (57) 

The parameter set Ja [Eq. (18)], with three-sublattice 
order, corresponds to p = 3. 

Our preliminary investigation of this question, using 
classical Monte Carlo simulation^^ and exact diagonal- 
isation^^, suggests that the answer is not simple. Eor 
p = 4, we do indeed find four-sublattice order. However 
thermal fluctuations favour a collinear state with “up- 
up-down-down” stripe order over the expected ID spi¬ 
ral. This collinear state cants in applied magnetic field, 
in a way reminiscent of the planar (XY) modeF^. Quan¬ 
tum fluctuations select the same, collinear, stripe order 
for p = 4 Ref. 45. Meanwhile, exact diagonalisation car¬ 
ried out for J = (—1, 0.3, 0.175), i.e. p = 6 , for clusters of 
up to Y = 36 spins, indicate that the ground state is a 
commensurate ID spiral. They also suggest the possibil¬ 
ity of a 1/3-magnetisation plateau in applied magnetic 
field. Clearly, this issue merits further investigation. 

B. Spin liquid and “ring” correlations 

The phase diagram shown in Eig. 8 , calculated for the 
parameter set J^, where different spiral phases meet, 
looks superficially similar to the one for the parameter 
set associated with three-sublattice order J^, shown in 
Eig. 2. However, the addition of a degenerate manifold 
of classical ground states leads to a number of crucial 
differences, especially at high temperature, and for low 
values of magnetic field. 

The most obvious consequence of the enlarged ground- 
state manifold is the “ring” seen in the structure factor, 
as shown in Eig. 9. This ring is present in both the low- 
temperature ordered phases, and the high-temperature 
paramagnet, for all values of magnetic field. Simulations 
for more generic parameter sets^^ confirm that this “ring” 
is a generic feature of the classical phase boundary be¬ 
tween ID and 2D spiral ground states [cf. Eig. 4]. 

At finite temperature, the number of low-lying states 
contributing to the ring in the structure factor, T^ji 
[Eq. (44)], scales with the size of the system [Eig. 10, 11]. 
As a consequence, the ring makes a macroscopic contri¬ 
bution to the entropy of the system, even though the 
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number of classical ground states, defined by Eq. (20), 
only scales with the linear dimension of the system L. 
The high-temperature paramagnet phase in Fig. 8 can 
therefore be viewed as a spin liquid with a finite mag¬ 
netisation. 

The spin liquid found in this two-dimensional frus¬ 
trated ferromagnet differs from well-known three- 
dimensional examples, such as spin ice^^, or the py- 
rochlore Heisenberg antiferromagnet^^ in that the classi¬ 
cal ground-state manifold is sub-extensive. It has more 
in common with problems where ground-state manifold 
correspond to lines or surfaces in reciprocal space, such 
as the chiral spin liquid found on the diamond lattice^^, 
or the “ring” and “pancake” liquids reported in the frus¬ 
trated honeycomb antiferromagneb^^. We note, however, 
that unlike the honeycomb lattice example, we do not 
find any crossover to a conventional paramagnet at high 
temperatures. 

In the absence of a more complete description, we refer 
to the new state as a “helicoidal” spin liquid, in recog¬ 
nition of the fact that the system displays short-range 
helicoidal correlations^^ — cf. Fig. 10. The better char¬ 
acterisation of this phase remains an open problem. 


C. Conical Spiral 

Besides supporting a spin liquid, the ground state man¬ 
ifold for the parameter set J b [Eq. (22)] has the special 
property that it contains a large number of simple, com¬ 
mensurate wave vectors, Qr [Eq. (42)]. This sets the 
stage for a rather unusual form of order-by-disorder , 
in which both single-Q and multiple-Q ordered states 
are selected by fluctuations from states belonging to the 
“ring” of allowed Q [cf. Fig. 9]. Both of these phases 
are present for h = 0, and persist for a finite range of 
magnetic field, as shown in Fig. 8. 

At the lowest temperatures, thermal fluctuations se¬ 
lect an incommensurate single-Q state with the charac¬ 
ter of conical spiral, familiar from studies of the Heisen¬ 
berg antiferromagnet on an anisotropic triangular lat- 
tice3,50,74,75 spiral has uniform magnetisation 5'^, 

and wave vector very close to Qig [Eq. (48)]. In the limit 
T ^ 0, we can use a low-temperature expansion to show 
that this incommensurate spiral makes a larger contri¬ 
bution to entropy than any other state in the ground- 
state manifold of the model (cf. Appendix B). The same 
state is found in classical Monte Carlo simulations at low 
temperatures T < O.ITJa and low values of magnetic 
field h ^ Ja — cf. Fig. 8. At low temperatures, but 
larger values of magnetic field, this incommensurate non- 
coplanar state gives way to the commensurate, coplanar 
“Y—state” familiar from the triangular-lattice antiferro¬ 
magnet (cf. Fig. 8 with Fig. I of Ref. 62). 

It is a widely quoted “rule of thumb” that fluctua¬ 
tions in frustrated magnets, whether quantum or clas¬ 
sical, favour coplanar — and if possible, collinear — 
states^^. It is also frequently assumed (without proof) 


that quantum and thermal fluctuations will favour the 
same state. The present case provides an interesting 
counter-example. Somewhat surprisingly, thermal fluc¬ 
tuations select an incommensurate spiral for T ^ 0, as 
discussed above. Meanwhile, as shown in Appendix D, 
quantum fluctuations, treated at the level of linear spin- 
wave theory, prefer a state with commensurate wave vec- 

tor QL^ub [Eq- (11)]- 

By analogy with the quantum Heisenberg antiferro¬ 
magnet on a triangular lattice^^, we anticipate that 
the resulting quantum ground state will support three- 
sublattice, coplanar, “120-degree” order for h = 0, giv¬ 
ing way to a coplanar “Y—state” in applied magnetic 
field. Thus, for low values of magnetic field, quantum 
and thermal effects compete, setting up the possibility 
of a phase transition from the commensurate quantum 
ground state into the incommensurate state favoured by 
thermal fluctuations, as temperature is increased. In 
practice this would mean that the coplanar “Y—state”, 
already favoured by thermal fluctuations for larger value 
of field — cf. Fig. 8 — would displace the conical spiral 
for T ^ 0. 


D. Vortex Crystal 

While the conical spiral is a relatively conventional 
phase, the structure of the vortex crystal marks it out 
as very different from anything found in triangular- 
lattice Heisenberg antiferromagnets. The vortex crystal 
is formed through the coherent superposition of classi¬ 
cal ground states with four, symmetry-related wave vec¬ 
tors [Fig. 14(b)]. It is not selected by harmonic thermal 
fluctuations [cf. Appendix C], but rather by anharmonic 
(interaction) effects which become important at higher 
temperatures. At a more intuitive level, this 4(5 state 
can be visualised as a lattice of vortices in plane, 

with alternating vector chirality [Fig. 1(c)]. The result¬ 
ing “crystal” has a lO-site unit cell [cf. Fig. 14(e)], and 
meets the usual definition of a magnetic supersolid, since 
it spontaneously breaks the translational symmetry of 
the lattice, while at the same time breaking spin-rotation 
symmetry about the direction of the magnetic field (cf. 
Refs. 7 and 62 and references therein). However it differs 
from other examples of magnetic super solids, in that it 
has a uniform magnetisation 5'^, and translational sym¬ 
metry is instead broken by the spin texture perpendicular 
to magnetic field. 

Crystals formed of vortices have been widely studied 
in the context of type-H superconductors, where they 
are known as a “vortex lattice”Lattices of magnetic 
vortices are, however, very unusual, recent work on Mott 
insulators notwithstanding^^. The 4(5 state, described 
above, is not an exact analogue of a superconducting 
vortex lattice, since it contains vortices with alternating 
circulation. None the less, it is an interesting step in that 
direction. The vortex crystal also proves to be a robust 
state — we have checked in simulation that it survives for 
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deviations in exchange parameters 6J ^ 0.01 |Ji|. This 
robustness reflects the fact that it is only possible to form 
a 4(5 state from classical ground states specific to the pa¬ 
rameter set Jb , ruling out any continuous transformation 
into an incommensurate state. 

Recently, more attention has been paid to a different 
kind of crystals of topological defects: lattice states of 
Skyrmions^^’^^’^^’^^. In a Skyrmion lattice the spin tex¬ 
ture also also winds in alternating directions, but the core 
of each Skyrmion is associated with a modulated . In 
fact, modulation along the direction seems to be a 
generic feature of multiple-Q states found in magnetic 
models^^’^^. This is a strong constraint, and might ac¬ 
count for the lack of other examples of vortex crystals in 
spin models. In the light of the recent work on Skyrmion 
lattices, it might also be interesting to explore how a vor¬ 
tex crystal couples to itinerant electrons. 

It would also be of value to investigate how the IQ and 
4(5 states found in classical simulations of evolve 

in the presence of quantum fluctuations. On general 
grounds, we anticipate that the phases shown in Fig. 8 
should survive for large (quantum) spin. In principle, this 
could be investigated by self consistent mean-field meth¬ 
ods. The picture for spin-I/2 is however harder to assess. 
Exact diagonalisation studies reveal a very rich ground 
state phase diagram in the absence of magnetic field, and 
suggest a disordered ground state for the parameter set 
Jb [Refs. 45 and 46]. However, given the small cluster 
sizes available, and the large unit cell of the 4(5 state, it 
is impossible to rule out a vortex crystal ground state for 
spin S = 1 / 2 . 


VI. CONCLUSIONS 

In this Article we have explored some of the novel 
phases which arise in a simple model of a two- 
dimensional frustrated ferromagnet — the J 1 -J 2 -J 3 
Heisenberg model on a square-lattice, [Eq. (I)]. 

We have chosen to concentrate on the effect of thermal, 
rather than quantum fluctuations, using large-scale clas¬ 
sical Monte Carlo simulation, complemented by linear 
spin-wave theory, to determine the finite-temperature 
phase diagram of the model in applied magnetic field. 

Two distinct parameter sets were considered. The first 
of these, Ja [Eq. (18)], was chosen in order to favour 
three-sublattice stripe order. In this case we find a phase 
diagram, Eig. 2, remarkably similar to that of the Heisen¬ 
berg antiferromagnet on a triangular lattice^^. In applied 
magnetic field, a “ 120 -degree” ground state transforms 
first into a canted “Y-state” and then into a collinear 
I/3-magnetisation plateau, illustrated in Eig. 1(c), be¬ 
fore interpolating through a 2:1 canted phase to satura¬ 
tion. 

The second parameter set, Jb [Eq- (22)], was cho¬ 
sen to lie at the classical phase boundary between one¬ 
dimensional and two-dimensional spiral ground states. 
Here an enlarged ground state manifold manifests itself 


in a “ring” structure in the spin structure factor, Eig. 1(b), 
and in a much richer phase diagram, Eig. 8. In this case, 
the phases found include a high-temperature spin-liquid 
and a number of new ordered phases, as well as the “Y- 
state”, I/3-magnetisation plateau and 2:1 canted phase 
familiar from the Heisenberg antiferromagnet on a trian¬ 
gular lattice. In particular, for low values of magnetic 
field, a multiple-Q state with the character of a vortex 
crystal [Eig. 1(c)] competes with an incommensurate con¬ 
ical spiral state. These results are quite striking, and the 
subtle play of order by disorder which drives the finite- 
temparature phase transition between the incommensu¬ 
rate conical spiral and vortex crystal states would have 
been very difficult to anticipate a priori. 

The abundance of exotic phases found in this simple 
model suggests that there is much still to learn about 
frustrated ferromagnets. In particular, questions of com- 
mensurability, and of the role of quantum fluctuations, 
deserve further investigation. With respect to commen- 
surability, an obvious open question is whether vortex 
crystals, like the one studied in this Article, can form dif¬ 
ferent lattices, or lattices with different lattice constants. 
And quantum effects are expected to lead to an even 
richer phenomenology, stabilising new forms of magnetic 
order, as well as valence bond solid and spin-nematic 
states, with exact diagonalisation hinting at the possi¬ 
bility of a spin-liquid ground state^^. And, while the 
square-lattice Heisenberg model 1-L^^ provides a simple 
context for these questions, we anticipate that the same 
effects may be observed in a much wider range of systems. 

We conclude with a few comments on experiment. A 
number of magnetic insulators have been proposed as ex¬ 
amples of square-lattice frustrated ferromagnets. These 
include the vanadium phosphates Pb 2 V 0 (P 04)2 [16], 
BaCdV 0 (P 04)2 [17], and SrZnV 0 (P 04)2 [18], and 
the copper-based topotactic ion-exchange systems 
(CuCl)LaNb 207 [19], (CuBr)Sr 2 Nb 30 io, [20], and 
(CuBr)Sr 2 Nb 30 io [21]. The properties of the vana¬ 
dium phosphates in applied magnetic field appear to 
well-described by a simple J 1 -J 2 model^^’^^. However 
the much richer phenomenology of the cuprates suggests 
more complicated exchange interactions, which need not 
have the full symmetry of the square lattice^^. 

In the absence of a clear experimental validation 
from, e.g., inelastic neutron scattering, the Hamiltonian 
[Eq. (I)] is probably best regarded as a toy model 
capturing some of the the interesting features of quasi 
two-dimensional frustrated ferromagnets, rather than a 
complete description of the exchange interactions in these 
layered cuprates. None the less, this model does score two 
notable successes, in providing a route to both the 1/3- 
magnetisation plateau observed (CuBr)Sr 2 Nb 30 io [20], 
and the helical order found in (CuBr)Sr 2 Nb 30 io [21]. In 
the light of this, it would be very interesting to see further 
experiments on square-lattice frustrated ferromagnets in 
applied magnetic field. 
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Appendix A: Low-temperature entropy in the 
harmonic approximation and its relation to 
spin-wave theory 


In this Appendix we provide details of the spin-wave 
calculations used to explore the different ordered phases 
selected by fluctuations. We start by reviewing the low- 
temperature expansion of classical spin model, and show 
how it can be used to explicitly calculate the free energy, 
and thereby the entropy, associated with a given form of 
magnetic order, to leading order in T [Section A I]. 

We then show how equivalent results for the free energy 
can be obtained from the spin-wave dispersion obtained 
within linear spin-wave theory for the associated quan¬ 
tum model, using the conical spiral state as a worked 
example [Section A 2]. 

Einally, we derive explicit expressions for the entropy 
associated with the different forms of magnetic order en¬ 
countered in this work : (i) conical spirals in applied 
magnetic field [Section A3], (ii) coplanar A-sublattice 
states with in-plane magnetic field, such as the Y state 
[Section A 4], (hi) canted conical A-sublattice states such 
as the 4Q state in applied field [Section A 5]. 


1. Entropy within a classical low—temperature 
expansion 


The low-temperature thermodynamics of a classical 
spin model of the form Eq. (I) can be calculated through 
an expansion of the Hamiltonian in small fluctuations 
about the ground state (T = 0) configuration. Eor this 
purpose, one may proceed as follows. Let us consider a 
magnetic ground state with K sublattices. Without loss 
of generality, it is possible to select a reference frame such 
that, in this ground state, the spin on the sublattice 
(z/ = I, • • • , A) in the mth unit cell at position is 


written as 


" CX 

^v>,m 


Sm{e^^rn) Sm{(f>^^rn) 

Sy 

u,m 

= 

COs(0,.,m) 

oz 

_u,m_ 


_sin(6»^,m) cos((/)^,m)_ 


(Al) 


We then choose for each site (z/, m) a local frame 
(e^,, e^, 6;^), rotated from the reference frame, such that 
ez coincides with the ground-state spin direction 
and expand the spin components in this frame as 


G — 


Viy', 

I - 


2 

u,m 


-\-y^ 


(A2) 


where and denote transverse components of a 
spin deviation. One has, in the harmonic approximation, 

n = n^^) (A3) 

where is the energy of the ground state and is 
quadratic in spin deviations. For N spins with Nc unit 
cells {Nc = N/K): 


(A4) 

where ^ is a 2A-dimensional vector of spin deviations, 

= [xiq, • • • , xk,Nc, yi,i, • " , yK,Nc] (A5) 

and J\4 a symmetric 2N x 2N matrix. By Eourier trans¬ 
forming the spin deviations 

x{y)u,m = (A6) 

VJVc ^ 

Eq. (A4) becomes 

^(2)^1 ^ e^A4(k)|k (A7) 

IcGMbz 

where 

Ck = , • • • • • • ,^K,k]- (A8) 

Here Mbz denotes the magnetic Brillouin zone for the 
A-sublattice structure, and AJ(k) is a 2K x 2K Hermi- 
tian matrix. Diagonalizing AI(k) by a unitary transfor¬ 
mation, one obtains the form 

1 

^zy,kCf^,kCf^, —k: (A9) 

v' kGMez 


where are the eigenvalues of AI(k). The free energy 
per site at low T is written as 


N 


N 


-TlnT-r^ 


+ 0{T‘^) , 


(AlO) 
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where 


with 


1 

-^ =-(ln/t^,k) =’ (^ 11 ) 

U=1 k 

is the entropy per spin associated with harmonic spin 
fluctuations^^. We note that the evaluation of (In/^^y^k) 
does not require the computation of the eigenvalues of 
Al(k) but only of its determinant |Al(k)|, since 

1 

(111 ^zy,k) = 111 ^zy,k 

k _ u 

= (A12) 

k 


2. Alternative derivation of entropy from linear 
spin—wave theory for a quantum model 


It is also possible to evaluate the determinant |Al(k)| 
[Eq. (A12)], which determines the low-temperature en¬ 
tropy of classical spins [Eq. (All)], starting from a linear 
spin-wave (LSW) theory for quantum spins. To show 
how this works, we first derive the large-^S, LSW expan¬ 
sion of [cf. Eq (1)1, and then take a classical limit, 

setting S = 1. 

We consider the same ground state configuration with 
K sublattices as discussed in Section A 1. After the rota¬ 
tion to the local frames and a Holstein-Primakoff trans¬ 
formation of spin operators into bosonic creation and an¬ 
nihilation operators al^ and a^y^k, 


■ y572(a-, 

-i^/s/2{a, 

S-al 





(A13) 


the harmonic Hamiltonian for spin S' is written in the 
form 


where 


^qu 


"Hi?? = 1 E [ak^(k)ak - Ak 


(A14) 


A _ \J 


I k • • • 5 lll,-k: • • • : HK,—kj 


(A15) 


M(k) denotes a 2K x 2K matrix, and Ak is a scalar func¬ 
tion which determines the zero-point energy associated 

with a given form of order. 

(“2^ 

We can diagonalise Hqu [Eq. (A14)] using a (para¬ 
unitary) Bogoliubov transformation^^ 

ak = Tkbk , (A16) 

such that 

r^Z-iTk = I-i, (A17) 


I-l = 


Ik 0 
0 -Ik ’ 


(A18) 


and di K X K identity matrix Ik^ brings into a di¬ 
agonal form 


with 


^^u^ = ^E[bfWbk-Ak 

k 


!^(k) = /_ir-V_iM(k)rk = 


cjk 0 
0 cj-k 


where c^k is the diagonal matrix 


CJk = 


^i,k 


^2,k 


<^K,k 


(A19) 


(A20) 


(A21) 


and cjjy^k is the spin wave frequency in the u branch at 
wave vector k. 

We now consider a classical limit in Eq. (A14) by set¬ 
ting S = 1, in order to compare it with the classical 
harmonic Hamiltonian Eq. (A4). In the classical limit 
the spin operators a^y^k and aj, ^ are replaced with com¬ 
plex conjugate scalar fields T/^jy^k and '0*_k- (This cor¬ 
responds to considering a path integral formulation and 
omitting the imaginary-time dependence in the fields.) 
In this classical limit, the Holstein-Primakoff transforma¬ 
tion Eq. (A13) of S' = 1 spins in the lowest order becomes 
equivalent to the expansion Eq. (A2) of classical spins, if 
we set 


V^zy,k — (^zy,k •) 

C,-k = ■ (A22) 

Thus the classical limit of the quantum S = 1 har¬ 
monic Hamiltonian becomes equivalent to the classical 
harmonic Hamiltonian under the transformation 


^k = T^k, 


with 


V’k = [^l.k’ • • • ’ ^KM^ ^l.k> • • • > V’if.k], 

and 
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Since Ak originates from commutation relations of opera¬ 
tors, it can be neglected in the classical limit. Comparing 
both classical Hamiltonians, one finds that 

i\M{k)ik = (A26) 

which entails 

M{k) = 'pM{k)T. (A27) 

From this relation one immediately obtains 

|A^(k)| = |M(k)|, (A28) 

since |T^||T| = 1, and from Eq. (A20) 

K 

|M(k)| = \n{k)\ = (A29) 

V 

Setting S' = 1, we obtain the relation 

K 

i>i(k)i=n-. ,k^zy, —k: (A30) 

V 

With this result in place, the entropy Sgw entering into 
the classical free energy T [Eq. (AlO)], can be calcu¬ 
lated from the linear, quantum spin-wave dispersion cj^y^k 
[Eq. (A21)], through 

^ = -(InKk) = ^ ^ln|7W(k)| 

k 

1 ^ 

=-—=-(Inwi^k) • (A31) 

k 

We will explicitly show in the next subsection that 
Eq. (A30) holds for spiral states. We note that the spin- 
wave frequencies could also be obtained by solving 
the semi-classical equations of motion for the quantum 
spin-1 model. 

In conclusion, the classical entropy per spin, Ssw/N^ 
can be calculated from either the determinant of the 
matrix Al(k) [Eq. (A12)] found within a classical low- 
temperature expansion, or the values of the dispersion 
<^zy,k obtained within LSW theory [Eq. (A21)]. The first 
approach is slightly simpler but the latter is also interest¬ 
ing as it provides information on the dynamics both for 
the classical and the quantum model in the semi-classical 
approximation. It may also be more rapidly implemented 
if the LSW theory of the model has been derived previ¬ 
ously. 


3. Conical spirals 

Eor the case of conical spirals, it is convenient to choose 
the reference frame such that the magnetic field h is along 
the axis, instead of the axis. The ground-state 
spins of the conical spirals have then the same projection 


on the Sy axis, whereas their projections on the S^-S^ 
plane describe a spiral with wave vector Q. So, they 
can be written as in Eq. (AI), with equal 0 for all sites, 
introducing just one sublattice {K = I) as = (prn = 
Q -R^, where is the location of the site m. After the 
rotation to the local frames (e^,, and the expansion 

in spin deviations, the harmonic Hamiltonian reads 


^(2)^ l^|Tk_^(k)|k 
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k) 2/—k] 
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M^^ik) My^lk) 
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Here A1(k) is the 2x2 Hermitian matrix with coefhcients 

M""(k) = l[J(k + Q) + J(k - Q)] - J(Q), 

Myy{k) = M^^{k) cos^ e + [J(k) - J(Q)] sin^ 6, 

M^yik) = (M^®(k))* = ^ [J(k + Q) - J(k - Q)] cos6», 

(A33) 


where J(k) is the Eourier transform of the interactions, 
given in Eq. (5). The canting angle 0 is fixed by 


cos 0 = 


h 

J{0) - J(Q) • 


(A34) 


Next we describe the LSW theory of the quantum 
spin-5' model. After the rotation to the local frames 
and a Holstein-Primakoff transformation of spin opera¬ 
tors into bosonic creation and destruction operators a^. 
and Uk, the harmonic Hamiltonian for spin 5 is written 


^qu = ^ E [4^(k)ak - Ak 

k 

5 

(A35) 

where 



^k “ (^k’ ^-k) 


(A36) 

and 



Ak = -5J(Q) 


(A37) 

[cf. Ref. 87]. The 2x2 matrix M(k) is given by 


. . _ rA(k) + C(k) B{k) 

H(k) A(k)-C(k) 

5 

(A38) 


with 

A(k) = |{(cos2 d + l)[J(k + Q) + J(k - Q)] - 4J(Q) 
+2 sin^ 6>J(k)}, 

B(k) = I sin^ 6»[J(k + Q) + J(k - Q) - 2J(k)], 

C(k) = I cos 6>[J(k + Q) - J(k - Q)] . (A39) 
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A Bogoliubov transformation 


ak = ^kbk 


brings 1-L 


( 2 ) 

qu 


into a diagonal form 


(A40) 


local frames, one obtains a harmonic Hamiltonian of the 
form Eq. (A4) with 


M(k) 


M^y(ky 


(A49) 


=^E[^kWbk-Ak 

k 


(A41) 


where M'^^ik) ^ = x^y) are K x K matrices given 

by 


with 


n{k) = 


cjk 0 
0 cj_k 


where the spin wave dispersion is 


jy/ (k) — Ju,u' (k) cos (/)iy^iy' \li cos Ny , 

Mll,(k) = J,y{k) + [hcOS0, - K] , 

(A42) = M^^Ck) = 0, (A50) 

with = (j)j^ — (pjyf and 


wk = v'[A(k) + B(k)][A(k)-B(k)] + C(k). (A43) 

In the case of 5 = 1 we have the following relations 

Af®®(k) = A(k) + B(k), 

Af^^(k) = A(k) - B{k), 

M^y{k) = iC{k) . (A44) 

One can indeed directly see from Eqs. (A43) and (A44) 
that 


= y] ( 0 ) cos (pjyy. (A51) 

v' 

There are K = 3 sublattices in the Y state. The ex¬ 
pressions for the angles are similar to those found for 
the triangular antiferromagnet 


01 = TT, 


02 


—03 = arccos - 



(A52) 


wkw_k = A^(k) - B^{k) - C^(k) 

= Af*“(k)Af^2'(k) - \M^y{k)\^ 

= |A4(k)|, (A45) 

and hence 

^ = -(Inwk) . (A46) 

k 


4. Spin wave dispersion in coplanar A-sublattice 
states 


We now turn to the case of coplanar ground states such 
as the Y state with an in-plane field (or the coplanar AQ 
state in zero field) which require the introduction of K 
sublattices. For the sake of greater generality we rewrite 
the Hamiltonian as 




{u,m]u' ,m') v^m 

(A47) 

The Fourier transform of the exchange couplings is de¬ 
fined by 


7 / fki — \ J r\ f ,m 


(A48) 


Choosing the reference frame such that the ground state 
spins lie in the plane, with h along the axis, 

the ground-state spins can then be written as in Eq. (Al), 
with ^ = 7r/2 for all sites and independent of 

the unit cell. After expressing the spin deviations in the 


where Ja is given in Eq. (13), whereas the exchange cou¬ 
plings Eq. (A48) are 

Ju,u(k^ — cos(A^y) “h t /3 cos( 2 A^y), — 1 , 2 ,3} 

Jl,2(k) = J2,3(k) = J3,l(k) 

= + 726*'=“= cos(fcj,) + :|e-2*'=“= . (A53) 

For the AQ states the number of sublattices is TC = 10 
and we do not show the matrix of couplings Jyy{k). 


5. Spin wave dispersion in canted A-sublattice 
states 

Here we consider canted conical states with K sub lat¬ 
tices. As for the case of conical spirals, the reference 
frame is chosen such that the magnetic field is along the 
axis and the ground-state spins can then be written 
as in Eq. (Al). These canted conical states have equal 
projections of the spins along the field (equal values of 
0 for all sites) whereas the spin projections in the S^-S^ 
plane form a K sublattice configuration and include, as a 
special case, the AQ state where the spins uniformly cant 
in the direction of the magnetic field. 

The ground state has equal 0 for all sites and (j)u,m = 4^v 
independent of the unit cell. The matrix AA{k) has the 
same structure as in Eq. (A49) but with 

Mll.{k) = J,y{k) 

,(k) = Mll,{k) cos" e + [J,.yik) - d^yN^] sin" 9, 
(1^) = ['^vy (k) - Jvy (-k)] sin cos 6, (A54) 
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where the canting angle 0 is related to the applied field 
by 


cos 0 = 


hK 

T.VV' - <^0S(t)uy) ' 


(A55) 


In the case of a commensurate spiral state, = Q • 
Rz.,0, one can check that, as expected, [Eq. (A4)] can 
be rewritten in the form Eq. (A32) with the coefficients 
Eq. (A33) where 


= . (A56) 

V 


Eor the same reasons as above, Eq. (A31) holds and the 
entropy factor can be alternatively computed from the 
spin-wave frequencies. 


6. Spin—wave dispersion in spiral states 

In the main text we investigate the fingerprints of the 
spin waves of spiral states in the uniform paramagnetic 
phase. Eor this purpose, the spin-wave spectrum of the 
fully polarized state provides useful information about 
the low-energy excitations above the ring. The dispersion 
relation is obtained at the saturation field in a straight¬ 
forward manner as 

a;(k)isL = J(k) - J(Q), (A57) 

which has the lowest energy at the spiral wave vector 
Q. By imposing a low-energy cutoff, we can identify 
momenta close to the ring, which contribute to the low- 
energy excitations, and use them in Eq. (44) to calculate 
the total structure factor associated with the ring. 

It would be interesting to see how this dispersion is 
modified in the spiral ordered phase. Here we consider 
the coplanar spiral state at zero field. The dispersion is 
obtained from Eq. (A43) by setting S = 1 and 0 = nj2 

w(k) = yi[J(k + Q) + J(k - Q) - 2J(Q)] [J(k) - J(Q)]. 

(A58) 

When k is close to the spiral wave vector Q, one finds 

a;(k) « v'J[J(k) - J(Q)], (A59) 

where J = [J(2Q) + J(0)]/2 — J(Q). At the parameter 
set Jb, the spin-wave dispersion has gapless line modes 
along k = [Eq. (20)] in the LSW approximation. 



FIG. 15. (Color online) Entropy per spin Ssw/N of the two- 
dimensional spiral states which make up the classical ground state 
manifold of [cf. Eq (1)] for the parameter set Jb [Eq. (22)] 

The inset shows the set of wave vectors associated with the ground- 
state manifold, as defined by Eq. (20). Thermal fluctuations select 
the state with highest entropy, which corresponds to a an incom¬ 
mensurate spiral with wave vector close to (7r/2,7r/3). Ssw/N was 
calculated through Eq. (A46), using the linear spin-wave dispersion 
C(;(k) [Eq. (A43)], as discussed in Appendix A. 


at low temperatures and vanishing magnetic field, as sug¬ 
gested by Monte Carlo simulations. In Eig. 15 we plot 
the the entropy per spin Ssw/N [Eq. (All)] for the the 
family of conical states with wave-vector Q = {Qx^Qy) 
satisfying the ring equation Eq. (20), in the absence of 
an applied magnetic field. This was computed from the 
spin-wave frequencies [Eq. (A46) and Eq. (A43)], and is 
Ssw/N plotted as a function of Qx- 

Thermal fluctuations select, within the degenerate 
ring, the spiral state with the highest entropy, which 
is described by an incommensurate wave-vector slightly 
away from (7r/2,7r/3), or wave-vectors related to it by 
symmetry. This result is supported by the Monte-Carlo 
simulations [see Sec. (IV B)], and continues to hold for a 
small applied field, where the state evolves from a copla¬ 
nar into a conical spiral. 


Appendix B: Entropy selection of a spiral with 
Q (7r/2,7r/3) out of the degenerate ring in the 
harmonic approximation 

In this Appendix we show how, for the parameter set 
Jb [Eq. (22)], the harmonic analysis developed in Ap¬ 
pendix A predicts the entropic selection of a spiral state 


Appendix C: Entropy competiton between the 
Q (7r/2, 7r/3) IQ cone and the coplanar Y states in 
the harmonic approximation 

In this Appendix we show how, for the parameter set 
Jb [Eq. (22)], a canted Y state prevails over a conical 
version of the spiral state at low temperatures, for a suf¬ 
ficiently strong applied magnetic field. In Eig. 16, we plot 
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FIG. 16. (Color online) Entropy per spin Ssw/N [Eq. (A46)] 
of different competing states as a function of magnetic field h. 
The entropy of the coplanar Y state increases with increasing field 
and, for 0.8Ja < h < overtakes that of the conical state with 

Q = (7r/2,7r/3). Entropy was calculated with linear spin-wave 
(LSW) theory, for the parameter set Jb [Eq. (22)], as discussed 
in Appendix A. Within this linear, low-temperature approxima¬ 
tion, the 4Q state is not selected by fluctuations at any value of 
magnetic field. 



FIG. 17. (Color online) Quantum correction Ae/S to the ground- 
state energy (per spin) of spiral states [Eq. (D4)] for the spin-AS 
quantum model discussed in Appendix A 2. Calculations were car¬ 
ried out within the linear spin-wave theory (LSW), for coplanar 
sprial states with wave vector Q belonging to the classical ground- 
state manifold in the absence of magnetic field, for the parameter 
set Jb [Eq. (22)]. The inset shows the range of wave vectors con¬ 
sidered. Quantum fluctuations select a sprial with Q = (27r/3,0), 
and symmetry-related states. 


the entropy per spin Ssw/^ [Eq. (All)] for a selection of 
different states within the degenerate ground state mani¬ 
fold, as a function of magnetic field, calculated using the 
results of Appendix A3 and Appendix A5. The entropy 
factor of the coplanar Y state increases with increasing 
field, until it is finally selected by at low-temperature 
fluctuations for O.SJa ^ h < 3Ja as supported by the 
Monte-Carlo simulations [see Sec. (IV B)]. The 4Q state 
is not favoured by thermal fluctuations at low temper¬ 
ature, which is again in agreement with the simulation 
results. However, its entropy factor at low field is rather 
close to the one of the most favored conical spirals. This 
small difference is overcome by anharmonic effects which 
stabilize the 4Q state at larger T for field h < O.SJa- 


Appendix D: Quantum selection out of the 
degenerate ring in the quantum model in the 
large— aS limit. 

In this Appendix we discusses the role of quantum fluc¬ 
tuations in selecting an ordered ground state from the 
degenerate ground state manifold for the parameter set 
Jb [Eq. (22)], within the linear spin-wave (LSW) theory 
developed in Appendix A 2. From the LSW expansion 
of the quantum model [Eq. (A19)], one finds that the 


ground-state energy per spin in the large-AS limit can be 
written as 


^qu — ^0 + Ae, 


(Dl) 


where eq is the ground-state energy per spin of the clas¬ 
sical model and Ae is its first-order correction inl/S due 
to the the zero point motion of the spin waves 


k 


(D2) 


Here Ae is proportional to S whereas eq is proportional 
to S‘^. Among a set of classically degenerate ground 
state, quantum fluctuations will favor the state which 
minimizes Ae. 

For conical spirals 

C2 

eo = ^ [ J(0) cos^ 0 J{Q) sin^ O] + hS cos 0, (D3) 

which reduces to eq = 5'^J(Q)/2 in zero field, and 

^e=^^co^ + p{Q). (D4) 

k 

In Fig. 17 we show Ae/S at the Jb point in zero 
field for the family of conical states with wave-vector 
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Q = {Qx^ Qy) satisfying the ring equation Eq. (20). One 
sees that quantum fluctuations favor the three-sublattice 
stripe states with Q = (±27r/3,0) or Q = (0,±27r/3). 

A similar calculation based on the LSW theory for the 
iC-sublattice cases [Eq. (A19)] shows that, in zero field, 


the 4Q state has a higher energy than the Q = (27r/3, 0) 
spiral state (or equivalently the Y state at zero field). In 
a finite field, the Y state becomes favored over the conical 
spirals and the canted 4Q state. 
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